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We study a layer of grains atop a plate which oscillates sinusoidally in the direction of gravity, using 
three-dimensional, time-dependent numerical solutions of continuum equations to Navier-Stokes 
order as well as hard-sphere molecular dynamics simulations. For high accelerational amplitudes of 
the plate, the layer exhibits a steady-state "density inversion" in which a high-density portion of the 
layer is supported by a lower-density portion. At low accelerational amplitudes, the layer exhibits 
oscillatory time dependence that is strongly correlated to the motion of the plate. We show that 
continuum simulations yield results consistent with molecular dynamics results in both regimes. 
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Although experimental 0, and computational 0- 
evidence demonstrate the potential for hydrodynamic 
models to describe important aspects of granular flow, a 
general set of governing equations for granular media is 
not yet recognized 8, 9]. Several proposed rapid granu- 
lar flow models use binary, inelastic hard-sphere collision 
operators in kinetic theory to derive equations of motion 
for the continuum fields: number density n, velocity u, 
and granular temperature T [i"ol - fl2l |. As Eshuis, et al [1] 
stated in 2010, "The holy grail question in research on 
granular dynamics is [HI, [Tj] , To what extent can gran- 
ular flow be described by a continuum approach?" 

Density inversion, in which a low-density region near 
the bottom of a granular layer supports a higher-density 
region above it, has proven to be significant for the study 
of granular hydrodynamics. This phenomenon has been 
identified in vertically shaken layers @, ll5T - fl7l | as well 
as in layers flowing parallel to a surface, such as in 
gravity-driven flow down an incline [Lsl [l9T | . 

In their seminal investigation Q , Lan and Rosato stud- 
ied density inversion in vertically vibrated granular me- 
dia. A layer of grains with depth H and uniform diameter 
a atop a plate that oscillates sinusoidally with frequency 
/ and amplitude A will leave the plate at some time in the 
oscillation cycle if the maximum acceleration of the plate 
A (27r/) exceeds the acceleration of gravity g. 
The oscillating plate can be characterized by the dimen- 
sionless parameters T = a max /g and /* = fy^H/g. Lan 
and Rosato studied density inversion in such a system by 
conducting soft-sphere discrete element method (DEM) 
simulations and comparing these results to kinetic theory 
predictions of Richman and Martin [20[ . 

These continuum predictions did not account for the 
time dependence of the layer or the plate, but rather 
treated the oscillating plate as a source of thermal energy 
and assumed one-dimensional (1-D) steady-state density 
and temperature distributions as functions of height in 
the cell. To characterize the rate of kinetic energy input 
through shaking, they used the dimensionless RMS speed 
of the bottom plate VJ, = (27r/^4) / {^/2ag) as their con- 



trol parameter. It has since become common to instead 
use the dimensionless shaking strength [l5[ 
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In Fig. 5 of their manuscript, Lan and Rosato chose 
a fixed "dimensionless mass hold-up" m t = N ££ =2.5, 
where N is the number of particles in the layer, and C 
is the cross-sectional area of the bed. This corresponds 
to an average layer depth H w 4.3cr as poured. They 
investigated shaking strengths 5 = 8 and S = 50 by 
varying T and /* in direct proportion, while holding S 
constant. For each value of S, their DEM simulations 
exhibited two different regimes of behavior, depending 
on r (or proportionally, /*). 

In the high-r regime, DEM simulations exhibited a 
nearly steady-state density profile, in which a high- 
density region above the plate was supported by a high- 
temperature but low-density region near the plate. Their 
continuum model produced density and temperature pro- 
files consistent with this steady-state "density inversion" 
[7|. In the low-r regime, however, the motion of the 
layer in DEM simulations exhibited strong time depen- 
dence correlated to the motion of the plate. Similar time 
dependence at low acceleration is found in conjunction 
with shocks [2l[, patterns (22|, and compaction [23| of 
granular beds. Lan and Rosato found that their contin- 
uum model could not accurately describe these density 
profiles Q in this regime. 

Since that time, several studies have focused mainly 
on the high-r, steady-state regime, providing significant 
evidence that steady-state hydrodynamic solutions can 
accurately model this regime [I^| - [l7| . A further study 
investigated the transition from a time-independent Lei- 
denfrost state (in which a high-density crystalline cluster 
is supported by a low-density, fast-moving granular gas) 
to a convection state using linear stability analysis of 
time-dependent continuum equations [9j. 

We believe that the failure of kinetic theory to describe 
the low-r regime was primarily due to the assumption 
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of a steady-state solution. In our current manuscript, 
we conduct time-dependent continuum simulations and 
compare to molecular dynamics (MD) simulations. This 
comparison is analogous to the procedure used by Lan 
and Rosato, with two significant modifications. First, 
while Lan and Rosato compared 1-D, time-independent 
kinetic theory predictions to soft-sphere DEM simula- 
tions, we use a 3-D, time-dependent continuum simula- 
tion. We choose a continuum model [ll| that has previ- 
ously yielded results consistent with a hard-sphere MD 
simulation for shocks [2l[ and patterns [j| in granular 
layers, and compare to that same MD simulation. Sec- 
ondly, while Lan and Rosato modeled sinusoidal shaking 
along all three spatial axes, we model vertical oscillations 
commonly found in experiments 0, [l5| . We demonstrate 
that continuum simulations including time dependence 
are consistent with density profiles from MD simulations 
for both the high-F and low-r regimes. 

We numerically integrate a set of time-dependent con- 
tinuum equations to Navier- Stokes order in three spatial 
dimensions, as proposed by Jenkins and Richman (Tl| 
for a dense gas of frictionless (smooth), inelastic hard 
spheres. The granular fluid is contained between impen- 
etrable horizontal plates at the top and bottom of the 
container. The lower plate oscillates sinusoidally between 
height 2 = and z = 2A, and the ceiling is located at a 
height L z = 320<r above the lower plate. Periodic bound- 
ary conditions are used in the horizontal directions x and 
y to eliminate sidewall effects in a box with horizontal di- 
mensions L x = L y = 22a. We use the same grid sizes 
and numerical techniques as those used in p| [J. 

We compare continuum results to results from a fric- 
tionless, inelastic hard-sphere MD simulation. A ver- 
sion of this simulation including frictional interactions 
has been validated by comparison to experimentally ob- 
served standing wave patterns [22J, [24J , size segregation of 
binary mixtures [25| . and fluidization of granular mono- 
layers [26| . The horizontal dimensions of the cell are iden- 
tical to those in the continuum simulations, with periodic 
boundary conditions. The particles are constrained be- 
tween a bottom plate which oscillates sinusoidally be- 
tween z = and z = 2A, and a ceiling fixed at z = 320a. 

The layer depth as poured H = 4.30ct and the normal 
coefficient of restitution e = 0.90 are chosen to match 
those in Each simulation runs for 150 cycles of the 
plate to reach an oscillatory state; number density n 
and granular temperature T are then calculated at each 
height in the cell by averaging the following 100 cycles. 

We use values of F previously used in Q to investigate 
the time-averaged volume fraction v = ^na 3 as a func- 
tion of height for S — 8 and S — 50. In Fig.[TJ both values 
of S show maximum density near the plate for low accel- 
erational amplitudes T = 2.84 and T = 5.68, while higher 
values of V show density inversion. In all cases, signifi- 
cant volume fraction is found at at higher values of z for 
the higher shaking strength S = 50. Previous kinetic the- 
ory predictions were consistent with DEM simulations for 
r = 28.3 and T = 56.5, but displayed qualitatively dif- 
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FIG. 1. (Color) Time-averaged volume fraction v as a func- 
tion of height z/a (ordinate) for varying F corresponding to 
(a) 5 = 8, and (b) S = 50. Volume fraction is calculated at 
eight equal intervals during each oscillation cycle, and these 
eight times are averaged over 100 cycles of the plate. Data 
from MD simulations are represented by symbols, while con- 
tinuum data are shown as curves. Although the cell height is 
z — 320cr, this figure only displays a height up to z = 40er. 

ferent behavior than DEM simulations for T = 5.68 and 
T = 2.84 Q- By contrast, our time-dependent contin- 
uum simulations show volume fractions similar to those 
displayed in MD simulations for our entire range of T. 

We examine the time dependence of two sets of param- 
eters for S = 8 in Figd]a-h. One set (r = 56.5, /* = 6.6) 
corresponds to what Lan and Rosato called the "high 
T regime", while the other (r = 5.68, /* = 0.66) falls 
within the "low T regime" . 

For the low-r regime (Fig[5]a-d), the volume fraction 
and temperature profiles show significant oscillatory time 
dependence, correlated to the plate motion The high-r 
regime (Fig [5] e-h) shows nearly steady-state profiles, in 
which a higher-temperature but lower-density region near 
the plate supports a cooler, higher- density region away 
from the plate (z > 5a). Continuum simulations produce 
volume fraction and temperature profiles consistent with 
MD simulations in both regimes. 

Similarly, we examine v and T as functions of height 
for a low-r case (r = 5.68, /* = 0.26) and a high-r 
case (r = 56.5, /* = 2.6) corresponding to S = 50 
in Fig. [2] i-p. Interestingly, the temperature profile is 
strongly time-dependent in both cases. However, while 
the low-r case again shows significant oscillatory time 
dependence in volume fraction, the high-r volume frac- 
tion again shows a nearly steady state density-inverted 
profile. There is a larger discrepancy between tempera- 
ture predictions from continuum and MD for the low-r 
case in which the temperature is strongly anisotropic [3] ■ 
However, even in this case, the scalar temperature in the 
continuum simulation is comparable to that found in MD 
simulations, and volume fraction profiles are nearly iden- 
tical in both cases. 

These results indicate that time-dependent continuum 
simulations can capture a transition from an oscilla- 
tory, low-r regime to a steady-state density inversion 
at high r. This transition was absent in previous time- 



3 



S = 8, T = 5.68, /* = 0.66 S = 8, F = 56.5, /* = 6j 



S = 50, T = 5.68, /* = 0.26 5" = 50, F = 56.5, /* = 2.6 



41 

z/a 





(a) 


1 


ft=0.00 











(b) 

ft=0.25 




0.2 0.4 0.6 

v, T/(1 00 g a) 





(e) 




ft=0.00 




p 



(f) 

ft=0.25 





(g) 

ft=0.50 









(h) 




ft=0.75 




— continuum v 




—cont. 7/(1 OOga) 




• MDv 




» MD 7/(1 OOga) 


K7> 



0.2 0.4 0.6 

v, 7/(1 OOga) 



(i) 

ft=0.00 



(m) 

ft=0.00 



(k) 

ft=0.50 



(o) 

ft=0.50 





(I) 

ft=0.75 




— continuum v 
■■■cont. 77(1 OOga) 
• MDv 

» MD 7/(1 OOga) 





0.2 0.4 0.6 

v, 7/(1 OOga) 





1 ft=0.25 




* (n) 








L| ft=0.25 



























(P) 


\ 




71=0.75 






— continuum v 






---cont. 7/(1 OOga) 






• MDv 




» MD T/(100ga) 











0.2 0.4 0.6 

v, 7/(1 OOga) 



FIG. 2. (Color Online) Volume fraction v and dimensionless granular temperature T / (100g<r) as functions of dimensionless 
height z/a (ordinate) at four times ft in the oscillation cycle. The first two columns correspond to S — 8, with the left-most 
column (a-d) representing the "low-F regime" (r = 5.68, /* = 0.66) and the second column (e-h) representing the "high-r 
regime" (F — 56.5, /* = 6.6)0]. The final two columns correspond to S = 50, with the low-F regime (F — 5.68, /* = 0.26) in 
subplots (i-1) and the high-F regime (F = 56.5, /* = 2.6) in subplots (m-p). Curves indicate data from continuum simulations, 
points indicate MD data, and the oscillating bottom of the cell is shown as a solid horizontal black line. Due to the need for a 
sufficient number of particles for statistical calculations, T in MD simulations is only calculated when v > 0.01. 



independent continuum simulations Q. To examine the 
nature of this transition in our continuum simulations, 
we vary T for constant 5 = 8, and measure the height 
of the maximum volume fraction z max as a function of 
time over 100 cycles of the plate, sampled eight times 
per cycle. Taking a discrete-time Fourier transform of 
this time sequence allows us to calculate a power spec- 
trum as a function of response frequency f r within the 
range //TOO < f r < 4/, where / is the driving frequency. 
The dominant response frequency f r ,max is the frequency 
with the maximum power P ma x- We calculate the max- 
imum dimensionless power P^ ax — P max / 1 go and cor- 
responding dimensionless response frequency f* max = 
fr,max y/H/g, and plot them as functions of /* in Fig. [3] 
^max i s relatively high for low driving frequencies, indi- 



cating significant time dependence at f* max . However, 
Pmax tnen drops to near zero for /* > 2, indicating a 
transition to a nearly time-independent state. The dom- 
inant dimensionless response frequency lies along the line 
f* max = /* both below and above the transition, sug- 
gesting that the frequency response of the layer is iden- 
tical with the driving frequency through the transition. 
The transition occurs when the driving frequency exceeds 
a critical value, at which point, the amplitude of response 
drops significantly. 

To examine the effects of this transition on density- 
inversion, we vary T and plot z max during eight equally 
spaced times during a cycle in Fig. |H averaged over 100 
cycles. For both S = 8 and S = 50, low-T shaking pro- 
duces strong time dependence, as seen by a large stan- 



4 



CD 



D 

° 3 

3 

e 2 

1 



r.max 

A 10" 6 P* 



• • • • 



4 f* 5 

FIG. 3. Dimensionless maximum P^ OI of the power spectrum 
(circles) and the dimensionless dominant frequency fr, m ax at 
which this power is found (triangles) as functions of dimen- 
sionless driving frequency /* in continuum simulations. Pmax 
is reduced by a factor of 10 6 to fit on the same scale as f* tTnax ■ 
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FIG. 4. Dimensionless height of the maximum volume frac- 
tion Zmax/a in continuum simulations as a function of F (bot- 
tom axis) and /* (top axis), for (a) S = 8, and (b) S = 50. 
Zmax is calculated eight times during a cycle; points show the 
average of these times and error bars show standard deviation. 



dard deviation of z max . As V increases, the layer transi- 
tions to a steady state in which the height of the maxi- 
mum does not vary significantly with time. In both cases, 
the peak of the layer contacts the plate for low-r shak- 
ing, while the peak of the layer is far from the plate for 
high T. S = 50 yields consistently higher z max than does 
the lower shaking strength S = 8, indicating a more sig- 
nificant density inversion. For given V, the amplitude 
of the plate is much larger for S = 50 than for 5 = 8, 
which may account for the larger variations in height for 
S = 50. The transition also appears to occur at higher /* 
for S = 50, indicating that the critical frequency above 
which the layer becomes time-independent differs for dif- 
fering shaking strengths. 

We conclude that for a given shaking strength S, time- 
dependent continuum simulations exhibit a transition 
from a strongly time-dependent oscillatory state at low 
r to a density inverted state that is nearly decoupled 
from the oscillation of the plate at high V. While pre- 
vious studies only successfully used continuum theory to 
describe the high-r state, we find that time-dependent 
continuum simulations also describe the low-r state and 
produce results consistent with MD simulations in both 
regimes. The transition to a steady-state density inver- 
sion occurs at different driving frequencies for different 
shaking strengths; a future systematic study of varying 
shaking strengths could contribute to our understanding 
of this transition. 
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